AWS2021_22_web

load("C:/Users/wyatt/OneDrive - Colostate/MS Research/R_Git/CameronPeakFire/all_hr_AWS")

#import 15-minute data
fifteen_minute_uf_0804 <- read.table('C:/Users/wyatt/OneDrive - Colostate/MS Research/R_Git/CameronPeakFire/AWS/data/CR1000_7_Fifteen_Minute_0804.dat', 
  sep = ",", header = TRUE, skip = "1") %>%
  slice(., -(1:2)) %>%
  mutate(TIMESTAMP = ymd_hms(TIMESTAMP), TIMESTAMP = TIMESTAMP - hours(7)) %>%
  mutate_if(is.character,as.numeric) %>%
  mutate(site = 'uf') 


fifteen_minute_bf_0804 <- read.table('AWS/data/CR1000XSeries - new_Fifteen_Minute_0804.dat', 
  sep = ",", header = TRUE, skip = "1") %>%
  slice(., -(1:2)) %>%
  mutate(TIMESTAMP = as.POSIXct(TIMESTAMP, tz = "UTC"))# %>%  #, TIMESTAMP = TIMESTAMP - hours(7)) %>%
  mutate_if(is.character,as.numeric) %>%
  mutate(site = 'bf')
  
# create soilVUE DF
bf_0804_soil <- fifteen_minute_bf_0804 %>% 
  mutate(TIMESTAMP_LOCAL = as.POSIXct(format(TIMESTAMP, tz = "us/mountain", orgin = 'GMT', usetz = TRUE)),
         .after = "TIMESTAMP")%>% 
  filter(TIMESTAMP_LOCAL >= "2022-08-04 00:00:00") %>% 
  select(-c(3:30))

write.csv(bf_0804_soil, file = "C:/Users/wreis/OneDrive - Colostate/MS Research/Meetings/08092022/BF soilVUE/bf_0804_soil.csv")

bf_0804_soil_VWC <- bf_0804_soil %>% 
  select(c(1,2, starts_with("VWC"))) %>% 
  pivot_longer(., cols = c(starts_with("VWC")), names_to = "depth", names_prefix = "VWC_", 
               values_to = "VWC", values_transform = as.numeric) %>%
  mutate(depth = str_remove(depth, "cm_Avg"))

bf_0804_soil_Ka <- bf_0804_soil %>% 
  select(c(1,2, starts_with("Ka"))) %>% 
  pivot_longer(., cols = c(starts_with("Ka")), names_to = "depth", names_prefix = "Ka_",
               values_to = "Ka", values_transform = as.numeric) %>%
  mutate(depth = str_remove(depth, "cm_Avg"))

bf_0804_soil_T <- bf_0804_soil %>% 
  select(c(1,2, starts_with("T_"))) %>% 
  pivot_longer(., cols = starts_with("T_"), names_to = "depth", names_prefix = "T_",
               values_to = "T", values_transform = as.numeric) %>%
  mutate(depth = str_remove(depth, "cm_Avg"))

bf_0804_soil_BulkEC <- bf_0804_soil %>% 
  select(c(1,2, starts_with("BulkEC"))) %>% 
  pivot_longer(., cols = starts_with("BulkEC"), names_to = "depth", names_prefix = "BulkEC_",
               values_to = "BulkEC", values_transform = as.numeric) %>%
  mutate(depth = str_remove(depth, "cm_Avg"))

soil_data <- bf_0804_soil_VWC %>% 
  full_join(bf_0804_soil_Ka)%>% 
  full_join(bf_0804_soil_T)%>% 
  full_join(bf_0804_soil_BulkEC)

start_times <- which(soil_data$TIMESTAMP_LOCAL %in% as.POSIXct(c("2022-08-04 10:30:00", 
                                                                 "2022-08-04 12:15:00", 
                                                                 "2022-08-04 13:30:00")))

lims <- as.POSIXct(strptime(c("2022-08-04 09:00:00", "2022-08-04 15:00:00"), format = "%Y-%m-%d %H:%M:%S"))

# plot data
Ka_plot <- ggplot(soil_data, aes(x = TIMESTAMP_LOCAL, y = Ka, color = depth)) +
  geom_line() +
  geom_vline(xintercept = as.numeric(soil_data$TIMESTAMP_LOCAL[start_times]), size = 1,linetype="dotdash", color = "red") +
  labs(title = "Ka Timeseries", x = "Date") +
  scale_color_discrete(limits = c("5", "10", "20", "30", "40", "50")) + 
  scale_x_datetime(limits = lims, date_labels = "%H:%M", date_breaks = "1 hour")

Ka_plot

BulkEC_plot <- ggplot(soil_data, aes(x = TIMESTAMP_LOCAL, y = BulkEC, color = depth)) +
  geom_line() +
  geom_vline(xintercept = as.numeric(soil_data$TIMESTAMP_LOCAL[start_times]), size = 1,linetype="dotdash", color = "red") +
  labs(title = "BulkEC Timeseries", x = "Date") +
  scale_color_discrete(limits = c("5", "10", "20", "30", "40", "50")) + 
  scale_x_datetime(limits = lims, date_labels = "%H:%M", date_breaks = "1 hour")

BulkEC_plot

VWC_plot <- ggplot(soil_data, aes(x = TIMESTAMP_LOCAL, y = VWC, color = depth)) +
  geom_line() +
  geom_vline(xintercept = as.numeric(soil_data$TIMESTAMP_LOCAL[start_times]), size = 1,linetype="dotdash", color = "red") +
  labs(title = "VWC Timeseries", x = "Date") +
  scale_color_discrete(limits = c("5", "10", "20", "30", "40", "50")) + 
  scale_x_datetime(limits = lims, date_labels = "%H:%M", date_breaks = "1 hour")

VWC_plot

T_plot <- ggplot(soil_data, aes(x = TIMESTAMP_LOCAL, y = T, color = depth)) +
  geom_line() +
  geom_vline(xintercept = as.numeric(soil_data$TIMESTAMP_LOCAL[start_times]), size = 1,linetype="dotdash", color = "red") +
  labs(title = "T Timeseries", x = "Date") +
  scale_color_discrete(limits = c("5", "10", "20", "30", "40", "50")) + 
  scale_x_datetime(limits = lims, date_labels = "%H:%M", date_breaks = "1 hour")

T_plot

#BF 20 and 30 cm
#UF 50 and 60 cm

#SoilVUE Ka Data
all_soil_Ka <- all_hr %>%
  subset(site != "jw") %>% 
  select(c(TIMESTAMP, site,
             Ka_20cm_Avg, Ka_30cm_Avg, Ka_50cm_Avg, Ka_60cm_Avg)) %>%
  filter(TIMESTAMP > ymd_hms("2021-11-01 00:00:00")) %>%
  mutate(date = as_date(TIMESTAMP))

#Daily data
all_soil_Ka_daily <- all_soil_Ka %>%
  select(-TIMESTAMP) %>%
  group_by(date, site) %>%
  summarise(across(everything(), mean, na.rm = TRUE))

## `summarise()` has grouped output by 'date'. You can override using the
## `.groups` argument.

all_soil_Ka_daily_long <- all_soil_Ka_daily %>% 
  pivot_longer(., cols = c(starts_with("Ka")), names_to = "depth", names_prefix = "Ka_",
               values_to = "Ka", values_transform = as.numeric) %>%
  mutate(depth = str_remove(depth, "cm_Avg"))%>% 
  na.omit() %>%
  filter((site == "bf" & (depth == "20" | depth == "30")) |
           (site == "uf" & (depth == "50" | depth == "60")))

Ka_plot <- ggplot(all_soil_Ka_daily_long, aes(x = date, y = Ka, color = site, linetype= depth)) +
  geom_line() +
  labs(title = "Ka Timeseries", x = "Date", color = "(Site,", linetype = "Depth)")

Ka_plot

ggplotly(Ka_plot)

#SoilVUE T Data
all_soil_T <- all_hr %>%
  subset(site != "jw") %>% 
  select(c(TIMESTAMP, site,
             T_5cm_Avg, T_10cm_Avg, T_20cm_Avg, T_30cm_Avg, T_40cm_Avg, T_50cm_Avg, T_60cm_Avg, 
             T_75cm_Avg, T_100cm_Avg)) %>%
  filter(TIMESTAMP > ymd_hms("2021-11-01 00:00:00")) %>%
  mutate(date = as_date(TIMESTAMP))

all_soil_T_hr_long <- all_soil_T %>% 
  select(-date) %>%
  pivot_longer(., cols = c(starts_with("T_")), names_to = "depth", names_prefix = "T_",
               values_to = "T", values_transform = as.numeric) %>%
  mutate(depth = str_remove(depth, "cm_Avg")) %>% 
  na.omit()

#Daily data
all_soil_T_daily <- all_soil_T %>%
  select(-TIMESTAMP) %>%
  group_by(date, site) %>%
  summarise(across(everything(), list(mean = mean), na.rm = TRUE))

## `summarise()` has grouped output by 'date'. You can override using the
## `.groups` argument.

all_soil_T_daily_long <- all_soil_T_daily %>% 
  pivot_longer(., cols = c(starts_with("T")), names_to = "depth", names_prefix = "T_",
               values_to = "T", values_transform = as.numeric) %>%
  mutate(depth = str_remove(depth, "cm_Avg"))

T_plot <- ggplot(all_soil_T_hr_long, aes(x = TIMESTAMP, y = T, color = site, linetype= depth)) +
  geom_line() +
  labs(title = "T Timeseries", x = "Date")

T_plot

ggplotly(T_plot) %>% layout(hovermode = "x unified")

#BF 20 and 30 cm
#UF 50 and 60 cm

## `summarise()` has grouped output by 'date'. You can override using the
## `.groups` argument.

netRad <- ggplot(all_rad_daily) + 
  geom_line(aes(x = date, y = cumNR, color = site, linetype = "net")) +
  geom_line(aes(x = date, y = cumSWnet, color = site, linetype = "SWnet"))+
  geom_line(aes(x = date, y = cumLWnet, color = site, linetype = "LWnet"))

ggplotly(netRad)%>% layout(hovermode = "x unified")

netRad_med <- ggplot(all_rad_daily) + 
  geom_line(aes(x = date, y = cumNR_med, color = site, linetype = "net")) +
  geom_line(aes(x = date, y = cumSWnet_med, color = site, linetype = "SWnet"))+
  geom_line(aes(x = date, y = cumLWnet_med, color = site, linetype = "LWnet"))

ggplotly(netRad_med)%>% layout(hovermode = "x unified")

# Shortwave
sw_daily<- ggplot(all_rad_daily) +
  geom_line(aes(x = date, y= SWin_Avg_sum, color = site, linetype = "SWin")) +
  geom_line(aes(x = date, y= SWout_Avg_sum, color = site, linetype = "SWout")) +
  ggtitle("SW daily sum") +
  scale_color_manual(values = c("black", "green", "blue"))
ggplotly(sw_daily)

swin_daily_cum<- ggplot(all_rad_daily) +
  geom_line(aes(x = date, y= cumSWin, color = site)) +
  ggtitle("SWin cumulative daily sum") +
  scale_color_manual(values = c("black", "green", "blue", "yellow"))
ggplotly(swin_daily_cum)

swin_daily_ratio_filtered <- ggplot(rad_day_filtered) +
  geom_line(aes(x = date, y= SWin_ratio, color = site)) +
  ggtitle("SWin daily ratio (UF/BF)") +
  scale_color_manual(values = c("black", "green", "blue"))
ggplotly(swin_daily_ratio_filtered)

swin_ratio_boxplot<- ggplot(rad_day_filtered, aes(x = site, y= SWin_ratio, color = site)) +
  geom_boxplot() +
  ggtitle("SWin_ratio (UF/BF) daily boxplot - filtered") +
  scale_color_manual(values = c("green", "blue"))
ggplotly(swin_ratio_boxplot)

swin_daily_boxplot<- ggplot(all_rad_daily, aes(x = site, y= SWin_Avg_sum, color = site)) +
  geom_boxplot() +
  ggtitle("SWin daily boxplot") +
  scale_color_manual(values = c("black", "green", "blue"))
ggplotly(swin_daily_boxplot)

# Albedo
albedo_daily<- ggplot(all_rad_daily, aes(x = date, y= SWalbedo_Avg_median, color = site)) +
  geom_line() +
  ggtitle("Albedo daily median (between 1000 and 1400)")+
  scale_color_manual(values = c("black", "green", "blue"))
ggplotly(albedo_daily) %>% layout(hovermode = "x unified")

albedo_daily_ratio<- ggplot(all_rad_daily, aes(x = date, y= SWalbedo_ratio, color = site)) +
  geom_line() +
  ggtitle("Albedo daily ratio (UF/BF)")+
  scale_color_manual(values = c("black", "green"))
ggplotly(albedo_daily_ratio) %>% layout(hovermode = "x unified")

albedo_daily_boxplot<- ggplot(all_rad_daily, aes(x = site, y= SWalbedo_Avg_median, color = site)) +
  geom_boxplot() +
  ggtitle("Albedo daily boxplot")+
  scale_color_manual(values = c("black", "green"))
ggplotly(albedo_daily_boxplot)

albedo_ratio_boxplot<- ggplot(all_rad_daily, aes(x = site, y= SWalbedo_ratio, color = site)) +
  geom_boxplot() +
  ggtitle("Albedo ratio (UF/BF) boxplot")+
  scale_color_manual(values = c("black", "green"))
ggplotly(albedo_ratio_boxplot)

# longwave
sb = 5.67E-8

all_rad <- all_rad %>%
  mutate(lw_out_temp = sb*(AirTC_Avg+273.15)^4, diff = LWout_Avg - lw_out_temp, net_Avg = LWnet_Avg + SWnet_Avg)

lw_hr <- ggplot(all_rad) +
  geom_line(aes(x = TIMESTAMP, y = LWout_Avg, color = site, linetype = "measured")) +
  geom_line(aes(x = TIMESTAMP, y = lw_out_temp, color = site, linetype = "approx."))
ggplotly(lw_hr) %>% layout(hovermode = "x unified")



net_hr <- ggplot(all_rad) +
  geom_line(aes(x = TIMESTAMP, y = LWnet_Avg, color = site, linetype = "LW_net")) +
  geom_line(aes(x = TIMESTAMP, y = SWnet_Avg, color = site, linetype = "SW_net")) +
  geom_line(aes(x = TIMESTAMP, y = net_Avg, color = site, linetype = "net"))
ggplotly(net_hr) %>% layout(hovermode = "x unified")

net_daily <- ggplot(all_rad_daily) +
  geom_line(aes(x = date, y = net_Avg_sum, color = site))
ggplotly(net_daily) %>% layout(hovermode = "x unified")

# Permativity
all_perm <- all_hr %>%
  select(TIMESTAMP, site, starts_with("Ka")) %>%
  pivot_longer(., cols = c(3:11), names_to = "depth", values_to = "Ka")

perm <- ggplot(all_perm) + 
  geom_line(aes(x = TIMESTAMP, y = Ka, color = site, linetype = depth))

ggplotly(perm)

# Wind speed
wind_daily <- all_hr %>% 
  select(c(1,2,11,12)) %>%
  mutate(date = date(TIMESTAMP)) %>% 
  group_by(site, date) %>% 
  summarise(wind_mean = mean(WS_ms_Avg, na.rm = TRUE))

## `summarise()` has grouped output by 'site'. You can override using the
## `.groups` argument.

wind_speed <- ggplot(wind_daily, aes(x = date, y = wind_mean, color = site)) +
  geom_line()

ggplotly(wind_speed)

ggsave("wind_speed.png")

## Saving 7 x 5 in image

wind_speed_box <- ggplot(all_hr, aes(x = site, y= WS_ms_Avg, color = site)) +
  geom_boxplot()

ggplotly(wind_speed_box)

ggsave("wind_speed_box.png")

## Saving 7 x 5 in image